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We examine Bose-Einstein condensation (BEC) for particles trapped in a harmonic potential 
by considering it as a transition in the length of permutation cycles that arise from wave-function 
symmetry. This "loop-gas" approach was originally developed by Feynman in his path-integral study 
of BEC for an homogeneous gas in a box. For the harmonic oscillator potential it is possible to treat 
the ideal gas exactly so that one can easily see how standard approximations become more accurate 
in the thermodynamic limit (TDL) . One clearly sees that the condensate is made up of very long 
permutation loops whose length fluctuates ever more wildly as the number of particles increases. In 
the TDL, the WKB approximation, equivalent to the standard approach to BEC, becomes precise 
for the non-condensate; however, this approximation neglects completely the long cycles that make 

r^ ■ up the condensate. We examine the exact form for the density matrix for the system and show how 

O ' it describes the condensate and behaves in the TDL. 

S 

C3 '. I. INTRODUCTION 

-)— > 
C/2 . 

In recent years various groups have achieved Bose-Einstein condensation (BEC) in alkali gases. |l| The gases have 

been trapped magnetically in potentials accurately approximated by harmonic-oscillator wells. Because of these 

experiments there has been a deluge of theoretical papers on BEC in harmonic potentials, both for the ideal and 

' , ■ the interacting gases. However, even before the recent experimental achievements, there were a few theoretical 

^ [ considerations of bosons in harmonic potentials @~ ||] including one by the author in this journal. AJP has 

Q ■ recently pubhshed two papers on BEC in trapped gases. [g|' 1^ One of the pleasant aspects of this system is the 

CJ relative ease with which substantial theoretical advances can be made with this system, especially since the diluteness 

'~~' of the experimental gases often makes the ideal gas a fairly good approximation. The result has been some new 

I ' insights into the nature of BEC. 
^ : BEC occurs most fundamentally because of the symmetry of the wave function. This property is necessary because 

ly^ ^ bosons are identical particles, so that one cannot tell the difference between the arrangements, say, (1,2) and (2, 1). 
t^^ I This means that forms of the wave function that have the particles in differing permutation arrangements are equally 
'""' • likely. For example, for an ideal gas of just three particles in single-particle states (f>i might have, as one possible state, 
\0 ■ 

Q^ • ipabc{ri,r2,r3) = — [(I)airi)(l)b{r2)(t>c{r3) + (f>a{r2)(f>b{ri)<l)cir3) 

■^ ! +(t)a{ri)(t)b{r3)(t'cir2) + 4>a[r3)(j)b{r2)cj)c{ri) 

<^ ' + (t>a{r2)(t)b{r3)(t)c{ri) + 0a(r-3)0fc(ri)0c(r2)] (1) 
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Let us call the first term of this function the direct term; the next three are pair-exchange terms, with just two 
'^ ' particles interchanged; the last two involve triple exchanges with (1, 2, 3) — > (2, 3, 1) or (3, 1, 2). We might represent 
^ , this by a triangle with lines connecting the particles that permute and, because of this picture, call this a permutation 
O ■ "cycle" or "loop." (A pair exchange is a special case of a loop although it looks like a line.) Any permutation of a 
set of particles can be broken down into a combination of loops and singlets (particles that are not exchanged, e.g., 
particle 3 in the second term above). 

Feynman developed an approach to Bose systems when he adapted his path-integral view of quantum mechanics to 
statistical mechanics. [^ His work led to new understanding of superfluid helium and of the A- transition. [Q^ |l|] 
C^ ' The path-integral approach involves the use of the configuration representation, p(ri , . . . , tn ; r'j^ , . . . , rj,^ ) , of the density 
matrix as a central entity, a form often particularly useful in providing physical insight. For example, the diagonal 
element gives the probability of finding the N particles at positions ri, ..., tn. In the case of bosons, the symmetry of 
the wave function requires that the many-body density matrix include the effects of all possible permutations of the 
N particles. Since the density matrix is bilinear in the wave functions possible for the system, it includes not only 
direct terms (first term squared in our three-particle example wave function), but also cross terms among various 
permutation arrangements. Such cross terms become important when the thermal de Broglie wavelength is longer 
than the separation between particles so that particle identity becomes a serious issue; that is, there is a substantial 
overlap between an arrangement and a permutation of that arrangement. 



One might expect that small-scale (say, triple) exchange processes (as in the product of the first and fifth terms in 
the above wave function) could occur relatively easily when three particles chance to be near one another. What is not 
so clear is under what circumstances very large-scale (say, 50-particle) exchanges can occur as a unit. The occurrence 
of such grand "ring-around-the-rosy" exchange loops (the overlap of the direct and large-scale exchange terms in the 
wave function) is precisely what happens after the onset of BEG. This view is quite different from the usual one of an 
avalanche of particles falling into the lowest state - although that view is also valid. At BEG there is an "explosion" 
in the size of possible permutation-cycle events. Because their possibility is so clearly a result of particle identity - 
more evidently than the avalanche into the same state - an analysis of this point of view is a fruitful undertaking for 
the new intuition it provides. Finding this "loop-gas" relation between BEC and the growth of possible permutation 
cycles is the focus of this paper. 

The equilibrium properties of a gas, including the onset of BEG, can be determined from the partition function, 
which is the sum or integral over the diagonal terms of the density matrix. Thus the partition function and quantities 
derived from it can show the effects of permutation loops. Although our result that BEG is accompanied by a large 
increase in the possible size of permutation loops is itself not so new; what we show is that considering the partition 
function for particles in a harmonic-oscillator potential allows a particularly clear demonstration of the interpretation 
in terms of permutation loops. 

In his book, Feynman [ |lO| gives two derivations of the partition function of the homogeneous boson ideal gas in a 
box with periodic boundary conditions. The first is standard (p. 32 of his text); the second, arising from his interest 
in the density matrix approach, involves counting all possible permutation loops to compute the density matrix and 
then taking the trace to get the partition function (p. 62). (This approach was developed earlier by Matsubara, |l^] 
although this article is not referenced by Feynman.) The latter method accurately gives the partition function above 
the transition, but fails, in the homogeneous gas treated, to give all the insight for which one might hope, because an 
approximation that must be made to do the needed energy sums neglects the effects of the very long permutations. 
On the other hand, our exact treatment of the ideal gas in a harmonic potential is able to avoid this unfortunate 
approximation, include the condensate permutations, and thereby provide a rather useful bit of missing pedagogy 
concerning BEG. 

Because we present an explicit computation of the density of the gas (from our calculation of the density matrix), 
we are able to demonstrate a curious feature of BEG in a trap. The density is, of course, largest at the center of 
the trap. As the temperature is lowered, BEC finally occurs when the density at the center reaches the same critical 
value that it has for a homogeneous gas with periodic boundary conditions. From the loop-gas point of view, one sees 
that, at this density and temperature, relatively short permutation loops suddenly join to make large loops involving 
macroscopic numbers of particles. In a box, these loops are spread out uniformly. However, in the trap, the density is 
critical near the center, so the macroscopic loops occur in that region. Of course, that feature corresponds precisely 
to the fact that the ground state of the oscillator, into which condensation is taking place, is more localized than any 
other state. 

The density matrix itself provides several useful criteria for the occurrence of BEG, even when interactions are 
present. Gertain integrals over the density matrix must be macroscopic (of order N), for example. |13| For a homo- 
geneous fluid of particles in a box with periodic boundary conditions, the one-body reduced matrix, pi(r, r'), which 
is the A^-body matrix traced over all but the variables corresponding to a single particle, has the property |13| 

£»i(r,r') ^ no/V for |r-r'| ^oo (2) 

where uq is the condensate number and V is the volume. This property has been useful in Monte Carlo simulations 
mM of boson systems, e.g., superfluid ^He, for finding uq theoretically. However, if the particles are in a trap, such as 
a harmonic well, then gi{r, r') will approach zero, simply because of the nature of the potential, when |r — r'| — > oo, 
and other criteria must then be found. For actual experiments in a trap, the condensate is recognized by its sudden 
appearance as a compact cloud of particles that forms at the center of the trap. 

By looking in detail at what is going on in configuration space and at permutation cycles, rather than by concen- 
trating on the usual view of BEG as a collapse of many particles into the lowest state, one is able to arrive at quite 
different physical insights. Indeed such insights, first developed by Matsubara and Feynman, have had important 
applications in current research. 

In recent years Ceperley and co-workers have developed a powerful numerical procedure, based on Feynman's path 
integral approach, known as Path-Integral Monte Carlo (PIMG). |15LThc technique works well for dense systems such 
as superfluid helium, ||l^ and for interacting dilute gases as welL[ 16| ' p7| ^ [|l8[ The one-body density matrix limit 
above is no help for the trapped gases so that in recent PIMG simulations |16[^ [H no has been inferred by looking 
for the number of particles involved in long permutation cycles. These turn out, not coincidentally, to be localized 
near the center of the trap and, it is argued, constitute the condensate. In the ideal-gas analysis of this paper, such 
a useful connection between long permutation cycles and the localized condensate will arise in a simple analytic way 



without the complications of the numerical simulations needed when interactions are present. (Although PIMC is 
somewhat related to the present paper, it is not a focus of this paper and its approach will not be investigated. For 
more information on it we refer the reader to the Ceperley review. |l5[) 

In another recent related treatment [^ of the dilute Bose, all orders of permutation cycles are summed analytically, 
but while also including the effects of interactions in low order. 

II. STANDARD TREATMENT 

We first consider a standard textbook-type treatment of BEC in a harmonic well - even though no textbook actually 
covers this at present. The harmonic potential we use is written 



Uir) ^ -Uo [-) (3) 

where r is the radial distance from the center of the potential, Uo is the strength of the potential, and we have included 
a scale factor R that will become useful later. The energy levels are given in Cartesian coordinates as 

3 

Em^^my,ni, = huj{mx + my + m^ + -) (4) 

with nix, "niy, rriz — 0, 1, 2, 3, ... The angular frequency u in this equation is related to the potential parameters by 



where m is the particle mass. 

We find the properties of the BEC by writing the average particle number in the grand canonical ensemble 

L-^ /3[?iaj(m^+m„+r?Xj + |)-;^] _ i ^ ' 

711^ ,niy ^niz 

where (3 is the inverse temperature, and fx the chemical potential. 

The potential is macroscopic in size - all N particles fit in it, of course - but to take the thermodynamic limit 
(TDL) we need to increase the number TV while at the same time keeping the average density constant. The way we 
do that here is to weaken the potential while increasing N . We now use the scale factor in Eq.(ra); that is, we assume 
that N/R^ ^ Nu;^ is to be kept constant while we increase N. pO] Further discussion of the TDL that comes to the 
same conclusion is considered in Refs. 4, 7, 21, and 22. The last reference justifies this feature by showing that it 
is necessary in order to keep the free energy extensive. An equivalent result arises in our treatment below. We can 
define a characteristic temperature Tq that remains constant in the TDL by writing 



1/3 



, C/o . Uo fdy^ kTo 



where d = N/R^, k is the Boltzmann constant. From Eq.d^), we see that Tg is given by 



k \ m 

From Eq.(|^), we see that the harmonic oscillator states are very closely spaced for large N so that one would expect 
that it would be a good approximation to replace the sums in Eq.(p|) by integrals. We write 

/■OO /"OO /"OO 1 

N = no + I dnix / dm^ / dnix - 
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where we have let u = (3hujm,x, etc. and 



a = -/3huj — Pfi. 



(10) 



As usual in Bose systems the summand in Eq. (O) does not vary uniformly near z — u + v + w^Q (it can peak very 
sharply there), and so the integral approximation miscounts ng, which becomes macroscopic. Thus we have had to 
include it as a special term in Eq.(0). We next expand the integrand in the usual way and integrate over the three 
variables. 
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is one of the Bose integrals. ||22 
If we use Eq. (^ , we find 



N = no + N 



yJ ^3(a). 



(13) 



From this equation we see that our procedure for taking the TDL has led us to a properly extensive expression for 
the non-condensed particles. Any relation between uj and N other than uj^ ~ 1/A^ would not lead to this result. For 
non-zero a , we can satisfy this equation with a negligible value of no- However, as T falls, the value of a needed to 
satisfy the equation approaches zero (becomes of order 1/A^) and F3 reaches its maximum value of C(3), where C{n) 
is a Riemann zeta-function. Bose-Einstein condensation takes place at this point and we have 



no 



N 



for T < Tc 



(14) 



where Tc = ToC(3)^^/'^ is the condensation temperature. In Fig. 1 we plot the condensate fraction no/N as a function 
of T as gotten from Eq.(|lJ), and also from the exact expression Eq.(|g) evaluated for finite N — 1000. |Q 

Eq.(O) contains a sum over I. Just what is the physical significance of that sum? For the above derivation it 
resulted only from a mathematical trick invoked to do an integral, but in fact it has a fundamental meaning in terms 
of permutation loops as we will show in Sec. IV. 

There is a considerably different looking derivation that has been used in the literature for this problem. ^' g] 
It involves the WKB, or semiclassical, approximation for the density n(p,r) of particles having momentum p and 
position r, which is 



n(p,r) 



h3 e/3[pV2m+U(r)-A<] _ 1 



(15) 



The density of particles at position r is the integral of this over p, which can be shown, by an expansion similar to 
that done above, to be 



^(^) = M 



rfp 



e/3[pV2m+U(r)-p] _ ]_ 



) ^ ^3/2 



V h^ 






(16) 



1=1 



If we integrate this over r with the quadratic potential of Eq.(g), we will find just the second term on the right-hand 
side of Eq.(13). That is, the WKB approximation does not include the condensate, which must be added on specially 
as in our previous derivation. 



We have given the standard approach to BEC for the trapped gas case. In the remaining sections we examine how 
we can get considerably more insight into BEC by looking at the process from the point of view of permutation loops; 
that is, we examine the so-called "loop gas." But first we consider an exact approach that avoids replacing sums by 
integrals. 

III. AN EXACT APPROACH 

When we replaced the sums by integrals above to get an answer that would be applicable in the thermodynamic 
limit, we lost the correct description of the condensate, which was added back in. In the case of the harmonically 
trapped gas, we can avoid this approximation by doing the sums exactly. In Eq.ffl) we make the expansion of the 
integrand as done in Eg. (|ll|) : 

oo 
AT _ \ ^ \ ^ -l[3h^{m^+my+m^)—la (\7\ 

Unlike the case of a free gas in a box where the energies go quadratically with the quantum numbers, the sums over 
mx,my, and ruz can be done exactly. The result is 

^ = E (l„'e-^;3n.)3 (18) 

Such expressions have been obtained in the literature previously. |2^ 

To retrieve the approximations of Sec. I, we need to note that us is of order 1/A^^/'' so that, for I smaller than iV^/'^, 
a very large number in the TDL, the argument of the exponential in the denominator is small. Expand in powers of 
it and keep only the first term of the expansion to find 



This is identically the same expression found in the last section for the non-condensed gas. The condensate has been 
neglected because a is of order 1/A^ for T < Tc as can be seen from the expression for no : 

oo 

no - = y e-"' ^ - (20) 

e" -1 ^^ a ^ ' 

1=1 

where the last approximation holds only below T^- The size of a means that the I sum is significant to values up 
to order N. Since the approximation we used to get Eq.(|l9|) was valid only for I < N^'^, it is no surprise that we 
have made an error in deriving Eq.(|l9|). We can see how this occurs in Fig. 2, which plots as a function of I, for 
finite N = 1000, the exact summand of Eq.(n8), the approximate summand of Eq.(n9), and the summand of uq of 
Eq.(20). (The last two summands do not quite add up to the exact total because of higher order terms neglected in 
Eq.(19). These would become smaller in the TDL.) The fact that the standard integration approximation neglects 
large / terms again draws attention to the question of the physical meaning of these terms. In the next section we 
show how the sum represents permutations cycles. 

IV. LOOP GAS APPROACH 



l3] and then later by Feynman in connection with his 
Dcen discussed by Elser [M , who seems to have invented 



The "loop gas" view was first introduced by Matsubara 
path-integral method in statistical mechanics jl^l ; it has also 
the terminololgy "loop gas." These references all deal with the homogeneous system with periodic box boundaries. 
The term "loops" is used here as shorthand for "permutation cycles" or "exchange cycles." The reader will find that 
the derivation given here is more complicated than that of the last two sections, but the motivation is to achieve new 
physical insight into the relation between BEC and the growth of permutation cycles. 

The boson iV-body density matrix involves only symmetric states and is [|lO[ 



ps(ri,--rN;r'i,...,r^) = — ^pD(rPi, ...rp„; r'^, ...,rN) (21) 

p 

where the sum is over all permutations and 

pi3(ri,...rN;r'i,...,ri^)= ^ e-''E'V^i(ri, ...,rNM(r'i, ...rj^) (22) 

all states 

is the density matrix for distinguishable particles. The partition function is the trace of this matrix or 

^=JnT. /rfri...drN(rp,,...,rp,|e-'3"|ri,...,rN) (23) 

• p J 

where the variable rpj represents the coordinate of the particle interchanged with particle j in permutation P. The 
sum over permutations here is just the result of the possible cross terms between various permuted elements of the 
symmetrized wave functions caused by particle identity. By using the configuration representation of the density 
matrix, we are able to track how the symmctrization affects the results. 

Any iV-particle permutation breaks up into a number of smaller loop permutations. For example, for N — 7 we 
might have a 4-particle loop, a pair permutation (2-particle loop), and a single particle. In the 4-particle loop, we 
might have particle #1 interchanged with #2, #2 interchanged with #3, #3 interchanged with #4, and #4 with #1. 
Given a set of configurations containing all C(gi, g2, ■••) ways of having qi loops of length 1, q2 loops of length 2, etc., 
we can write 

^^ = M E C{q^,q2,...)Y[Z{ir (24) 

{91, «2---} l 

where the sum is over all combinations of permutation numbers such that 

J2qil = N, (25) 

I 

and 

Zil)= fdr^...dri{rp^,...,rp^\e-^^\ri,...,ri) (26) 

corresponds to a unbroken permutation loop of / particles. Feynman ||l^ has given an argument to show that this 
quantity is 

While the argument given by Feynman can be consulted, for completeness we rederive this result in the Appendix. 
Eq.( p4[), with the result Eq.(^7|), is quoted in a beautiful little paper on Bose and Fermi partition functions by D. I. 
Ford ]27||, as well as by Matsubara and Feynman. 

To evaluate the most probable distribution of loops, we invoke the grand canonical ensemble and multiply Zn by 
g/3^Ar ^^^ g^j_|^ Q.^gj. ^Y[ particle N values. This trick removes the restriction of Eq. (p5|) and allows the individual sums 
on qi to be done. The grand canonical partition function is 



2 



Z{l)e 



PtJ.1 \ *?' 



N I qi 
JJcxp^^^ (28) 



; 



What we want from this is the average particle number, (N) , and the average number of particles, (qi) , in the 
permutation loop of length I. We use the standard relation 

(7V)-|-fcTlnZ = ^Z(|)l/''^I (29) 



From this it is clear that the average value of the loop number qi is 

ill) - ^^ (30) 

Because H is the Hamiltonian corresponding to an ideal gas, Z{1) can be broken up into a product of single-particle 
densities. Z{1) is most easily evaluated by making each single-particle operator e~'^^' diagonal by interposing complete 
sets of harmonic oscillator functions as follows: 



Z{1)= Y^ / dri...dri(rp,,...,rpjmi,...,i 



.mi; 

mi ,7712,... ' 

xe-'^(^"-+-+^"<)(mi,...,m,|ri,...,ri) (31) 

where we have included only one subscript rui to represent all three quantum numbers for a single-particle state. The 
r integrations lead to factors like 

dri (mp^|ri) (ri|mi) = (5mp^,mi (32) 

where Pk represents the particle number that moved into position ri in the permutation. Since we have a ring 
exchange these (5-functions mean that all rrii become equal to one of them, call it m. The result is simply 



Z(/)=^e-^'^" (33) 



which is a one-body partition function at inverse temperature (31. The sum, over the three quantum numbers repre- 
sented by m is easily done: 



-?,pinuj/2 



Z(l) = T (34) 

The numerator comes from the zero point energy in each Em- Putting this result back into Eq.(|29|) (and dropping 
the averaging brackets), we find 



^ = E 



^ (35) 

7^ (1 - e-Pi^-'f 



which, as expected, is precisely the same result as derived in Sec. III. However, this derivation has led to the very 
instructive result that the sum over I is a sum over permutation loops, which arise from the symmetrization of the wave 
function. Note also that in the homogenous gas in a box, one is unable to do the exact sum over states as we could in 
Eq. (p3|) . In that case |QQ|^ [n2|^ [E6l one replaces the sum by an integral; the result is that the expression equivalent 
to Eq.(|35[) is approximate and neglects the condensate, which must once more be added in by hand. Eq.(^5|), on the 
other hand, is exact and includes the condensate. 

Figure 2 illustrates the fact that the condensate is made up of large-scale permutation loops, a fact that is well- 
known, but has not been so well-illustrated previously. What has perhaps not been so clear previously is that loops 
of all sizes contribute nearly equally to the condensate. The average number of particles pi — (qi) I in a permutation 



loop is given by Eq.(30) as 

Q-al 
P'- ^ n _ g-iTo/TAfi/3-j3 (36) 

where we have used Eq.(n). Because a is of order 1/A^ at T < Tc, the I sum extends to order N, but we know that 
the non-condensed particles are given by terms with ITq/TN^^^ <C 1, that is, for loops of size less than of order N^^^. 
The terms for larger I 3> N^'^ have pi = e~"', which is identical with the form given in the second last version of 
Eq.(po[) for riQ. All these terms then contribute only to the condensate. Each such term in the condensate is crudely 
of order 1, but there are of order N such terms so they contribute a total number of particles on the order of N. The 
fact that permutation loops of length of order I > N^^^ have only 1 particle in them on average means that the long 
permutation loops are forming and breaking with fluctuations in size that are extremely wild. |2q] 



V. DENSITY MATRIX 

The density matrix is a very useful quantity allowing averages of any thermodynamic variable to be evaluated. The 
one-body reduced density, mentioned in Sec. I, is particularly useful in considerations of BEC as Penrose and Onsager 



L3[ have shown. This quantity is defined by integrating the A'^-body density, Eq. (21) over all but one coordinate: 



£ii(r,r') = / dr2...drN ps(r,r2---rN; r', r2---rN) (37) 

The diagonal element Qi{r, r) is just the density of particles g{r) at position r. It turns out, after a short derivation, 
that the one-body reduced density for the boson ideal gas has the rather simple, physically intuitive form 

ei(r,r') = 5]n™^„(r)V';,(r') (38) 

m 

where 

^ (39) 



3/3(B„-m) _ 1 



is simply the occupation of the rnth single-particle state as used in Eq.(|6|), and ^m is the corresponding single-particle 
state. 

The density matrix gi{r,r') can be considered to be an operator with eigenvalues. In fact the spectral form of 
Eq.(pq) shows that the eigenfunctions are the V'm and the eigenvalues are the fim', it is straightforward to verify that 
/dr'gi(r,r')?/;in(r') = nm^m(i")- The fact that gi{r,r') has a large eigenvalue, in this case no, leads to its usefulness 
as a criterion for BEC. Even when there are interactions, the large eigenvalue is the condensate number and the 
eigenfunction is the condensate "wave function," but note that the solution is no longer simply the ground-state of 
the oscillator potential. In the case of weakly interacting bosons, this condensate wave function is approximately the 
solution of the well-known Gross-Pitaevskii nonlinear Schrodinger equation. |2| 

We can transform the density into a sum over permutation loops by expanding rim in the now familiar way. We get 

ei(r,r') = f;e^^'^'|^e-'^=='"^^(r)^l(r')| (40) 

1=1 I m J 

The quantity in curly brackets is simply the harmonic oscillator density matrix d{r, r') for a single free (distinguishable) 



particle, but at inverse temperature 1(3, rather than simply /3. The value of this is well known, and is |10 

3/2 



"'("■'^(iJI 



sinh {iphuj) ^ 



(41) 



X 

i—x^y.z 



One can see, for example, that d;(r, 0), and therefore Qi{r, 0), goes to zero at large r in contrast to the homogeneous 
case of Eq.(^). 

We are particularly interested in the particle density p(r) = gi(j, r). With a bit of hyperbolic function manipulation, 
one finds this to be 



One can integrate this p(r) to verify the result of Eq. (hi 

We have seen that the WKB approximation can be found by assuming that lf3fko is small. If we do that here we 
find 






This result agrees with Eq.dTo) for the harmonic potential. 



The largest value of the density occurs at the origin, of course, and we can ask what this value is at the BEC 
transition. For T above or equal to the transition temperature, the WKB approximation is valid in the TDL. Right 
at the transition the chemical potential vanishes (actually it is the quantity a of Eq.([l^) that vanishes, but the WKB 
approximation neglects the zero-point energy, which is of order 1/A^^/^.) We find that the critical density is 

MO)=(^^j E|^-(^^j C(3/2) (44) 

or 

*^— m U(3/2)j ^''^ 

But this is precisely the transition temperature for an ideal gas in a box at uniform density Pc(0)! Thus the trapped 
gas, which has its BEC localized near the origin, has a condensate there underjust the same conditions as the gas in 
a box. This interesting situation has been noted previously in the literature. ||]' Q' [P9| . 

In Fig. 3 we plot the full density as a function of position as well as the WKB approximation to it. One sees that 
the two are equal only for large r. In the case of A^ = 1000 at the temperature considered, the condensate number 
is about 456, so the non-condensed particle number is 544, despite the appearance in the plot that the condensate is 
much larger than the non-condensate. This feature is explained by the fact that there is a factor r^ involved when 
one integrates this to find the total N value; this factor diminishes the interior points and emphasizes the wide wings. 

Another way we can see that the high I permutation loops belong to the condensate is by considering the mean 
square radius of the individual permutation densities. This is 

2 /drrVi(r) 

CT; = "^3 rr- (46) 

/drpi(r) 

where pi is the summand of Eq. (B3) . The result is 

2 3 h 1 

^^ ~ 2 muj tanh I Phuj/ 2 ^ ^ 

The WKB approximation for this is gotten by considering l(3huj small, for which we find 

af(WKB) = 3 — -^ (48) 

This quantity drops off hke l/l. The contribution from the condensate comes from large iphcj {N^^^ < I < N) and is 

erf (condensate) — (49) 

2 rauj 

This is simply the mean square width of the ground-state of the harmonic oscillator. Figure 4 shows the plot of af , 
erf (WKB), and erf (condensate) for N — 1000. We again see that, for / > 40 or so, that the permutation loops all refer 
to the condensate. 

In recent path- integral Monte Carlo calculations [|l6[' [|8| for trapped gases, it was necessary to identify the con- 
densate by use of the permutation cycles. Ref.16 assumed that the size of the largest cycle occurring was equivalent 
to the condensate number. Eq] Ref.18, on the other hand, assumed that the condensate was made up of all particles 
for I greater than a value where af has become constant. However, we see that, for the latter approach, some of the 
particles having small permutation cycles also contribute to the condensate and must be counted. 

VI. DISCUSSION 

We have seen that it is possible and instructive to examine BEC via the permutation loop picture. Ideal gases 
trapped in a harmonic potential become particularly helpful in this regard because we can do many of the computa- 
tions analytically and exactly. The very long permutation loops are characteristic of BEC and form when the thermal 
de Broglie wavelength becomes comparable to the average distance between particles. Under these conditions multi- 
particle exchange becomes possible and the coherence characteristic of BEC takes place with the particles exhibiting 
the strong symptoms of indistinguishability. 



The best known characteristic of BEC is the presence of a macroscopic fraction of the particles in the lowest state. 
For the ideal gas, in a harmonic trap at absolute zero, those particles would be in the oscillator ground-state. From 
the point of view of permutations, this situation is equivalent to all possible lengths of permutation cycles being 
equally likely (up to length N, which is taken to infinity in the TDL). As the temperature is raised, but still below 
the condensation temperature, the shorter loops become increasingly prevalent. Shorter loops occur when particle 
indistinguishability is less a factor and corresponds to particles being spread out over excited oscillator states. Finally, 
above the transition temperature, no state has a macroscopic occupation of particles and most permutation cycles 
are short; very long ones are rare. 

BEC and the A-transition in liquid helium have been compared with a polymer transition, pq ] with long chain 
polymers forming below the transition temperature. In fact there is a transition in sulphur |3C| ] in which the opposite 
effect occurs: above the transition there are very long chains of sulphur atoms, while below the transition only eight- 
atom rings form; the isomorphism involves /3 —f kT. Interestingly the specific heat of sulphur has a "A-transition." 

The reduced single-particle density matrix has the lowest state as an eigenfunction with the condensate number as 
eigenvalue. For the ideal gas demonstrating this is a straightforward matter with its connection to long permutation 
loops made evident in the density. However, the property carries over into the interacting gas as well. Finding the 
eigenfunction of the one-particle density is one of the most important theoretical problems in the interacting fluid, 
with the solution representing the condensate wave function and the "order parameter" of the state. The solution 
of the Gross-Pitaevskii equation is a standard approximation to this wave function. In the case of real harmonically 
trapped gases, which actually do interact, this equation seems to represent the experimental situation fairly accurately. 

@ 
The PIMC approach was invented by Feynman and has been exploited in a wide range of fields, including the 

description of polymers. |3l| One of its most successful applications has been to bosons and liquid ''He by Ceperley 

and co-workers. [ [l5[ The approach considers in detail the particle positions in space, using the expression Eq.(|2l|) 

as its basis, and takes the formation of permutation loop "polymers" as one of its most important features. The 

interested reader should consult the review of Ceperley. |15( | 

APPENDIX 

In this Appendix we derive the expression Eq.(|27|) for the number of ways of forming qi permutation loops of 
length 1, (72 of length 2, ...qi of length /, etc. Suppose we consider an 9-particle system having two 3-particle loops 
[(1^ 2 ^ 3 ^ 1) and (4^ 5 ^ 6 -> 4)], a 2-particle loop [(7^ 8 ^ 7)], and a singlet [(9)]. We can mix up these 
particles to get other non-equivalent 3-, 2-, and 1-loops. We want to count ways we can do this. For example, suppose 
we replace 1 by 4 and 4 by 1, giving (4— > 2 ^ 3 — > 1), (1-^ 5 — > 6 — > 1), (7 — > 8 — > 7), and (9). This is a new 
arrangement. There are 5! different ways of mixing the particles to get new permutations arrangements, but not all 
of them are distinct. For example if we replace 2 by 1, 3 by 2, and 1 by 3, we get the permutation loop (3-^ 1^2^ 
3), which is just the same as the original arrangement, with the loop rotated one click, which makes no real change. 
In a loop of length I there are I such loop rotations that don't count for new arrangements. Another way of replacing 
particles that does not lead to new arrangements is by interchanging all the particles in one 3-loop of the same length 
by those of the other 3-loop. There are 2! such interchanges of loops possible in our special case. Thus the total 
number of possible distinct arrangements in our 9-particle case is 9!/3^2!. For the N-particle case we easily generalize 
to find that the number of distinct arrangements for the N particle case is 

where A'^! is the total number of ways of rearranging N particles, I'" is the number of ways of "rotating" the qi permu- 
tations loops of length I, none of which produces a new arrangement, and qi\ is the number of ways of interchanging 
whole loops of the same length, which also does not produce a new arrangement. The expression of Eq.(|50|) is known 
as "Cauchy's formula." [ p2[ 
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Figure 1. Condensate fraction versus T/Tq. The condensation transition temperature is To/C(3)^/'^ — 0.94To. The 
sohd hue is the thermodynamic hmit case given by Eq.(|lJ). The dotted hne is the exact solution for N =1000. 

Figure 2. Summands of various expressions of Sec. Ill for N versus summation index I. Here we have N ~ 1000 



and T/Tq = 0.7. The sohd hue is the summand of the exact expression Eq.(18); the dashed hne represents the WKB 
summand of Eq.(n3); and the dotted hne is the summand of the condensate in Eq.(pQ). The maximum value of the 
exact case at Z = 1 is 423. Since the index I corresponds to permutation loops, that means there are 423 singlets; the 
condensate (the dotted line) includes all the large I values (the long permutation loops) as well as some of the small 
loops and corresponds to a total of 456 particles; the remaining small, non-condensate loops contribute 121 particles. 
Note that the summands of Eq.([l9|) and Eq.(pCl|) do not quite add up to the exact result, because of the higher order 
neglected terms in the WKB expression. 

Figure 3. Density as a function of r. The solid line is the exact density, the dashed line is the WKB approximation, 
which represents the non-condensed particles, and the dotted line is the condensate density, derived using the first 
term in sum of Eq.(B8h. T/Tq = 0.8 here. Density is in units of [jnuj/fif/'^ and position in units of {Ti/mLoY^'^. 



Figure 4- Mean square width (in miits of (h/muj)) for particles in permutation loop / versus I. The solid line is the 
exact result, the dashed line is the WKB approximation to the mean-square width, representing the non-condensed 
particles, and the dotted line is the condensate width (3/2). 
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